function [f1,f2] = cal_Crack_RHS(E, nu, mode)
% calculate the RHS for Crack Problem using symbolic toolbox
% mode could be 1 or 2
%

syms r theta x y;


kappa = 3 - 4*nu;
G = E/(2 + 2*nu);

% for mode 1:
if mode == 1
 lambda = 0.5444837367825;  Q = 0.5430755788367;
 a = kappa - Q*(lambda + 1); 
 b = kappa + Q*(lambda + 1); 
 u = 0.5/G*((x*x+y*y)^(lambda/2))*(a*cos(lambda*atan(y/x)) - lambda*cos((lambda - 2)*atan(y/x)));
 v = 0.5/G*((x*x+y*y)^(lambda/2))*(b*sin(lambda*atan(y/x)) + lambda*sin((lambda - 2)*atan(y/x)));
else
    
% for mode 2:
 lambda = 0.9085291898461;  Q = -0.2189232362488;
 a = kappa - Q*(lambda + 1); 
 b = kappa + Q*(lambda + 1); 
 u =  0.5/G*((x*x+y*y)^(lambda/2))*(a*sin(lambda*atan(y/x)) - lambda*sin((lambda - 2)*atan(y/x)));
 v = -0.5/G*((x*x+y*y)^(lambda/2))*(b*cos(lambda*atan(y/x)) + lambda*cos((lambda - 2)*atan(y/x)));
end

fac = -E*(1 - nu*nu)/(1 - 2*nu);
f1 = fac*(diff(u,x,2) + (1 - 2*nu)/(2 - 2*nu)*diff(u,y,2) + 1/(2 - 2*nu)*diff(diff(v,x,1),y,1));
f2 = fac*(diff(v,y,2) + (1 - 2*nu)/(2 - 2*nu)*diff(v,x,2) + 1/(2 - 2*nu)*diff(diff(u,x,1),y,1));


%
%  [f1, f2] = cal_Crack_RHS(1, 0.3, 1);
% 